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Using the technique of Gaussian quadrature rules, a new estimator 
is proposed for approximating the distribution of a random variable 
given only a finite number of its moments. The estimator is shown by 
numerous examples to be accurate on the tails of both continuous and 
discrete distributions. Efficient algorithms exist for computing the 
estimator from the first 2N moments of the random variable. A robust 
implementation of the estimator is presented, along with rules that 
provide additional protection against computer roundoff errors. 

I. INTRODUCTION 

In this paper we present a method for computing the Cumulative 
Distribution Function (cdf) of an arbitrary random variable. Using 
the theory of Gaussian Quadrature Rules <gqrs), we derive an esti- 
mator that converges asymptotically to the true cdf. In practice, 
convergence is obtained without excessive computation. A general 
estimator is developed here that is applicable to a wide class of 
problems. 

Section 2.1 begins with a review of gqr analysis as it has traditionally 
been used for numerical integration. Several authors have shown the 
existence of extremely efficient algorithms for computing the param- 
eters of the gqr. An efficient and robust procedure for obtaining the 
gqr parameters is presented in the appendix. Two cdf estimators 
based on gqr are derived in Sections 2.3 and 2.4. The first estimator 
is most suited to numerical integration schemes and estimation of 
discrete distributions, while the second is appropriate for continuous 
distributions such as Gaussian noise or crosstalk. Section III gives 
numerous examples that show the inherent accuracy of the technique 
for continuous, discrete, and mixed distributions. Computational meth- 
ods for deriving the required moments are discussed, along with 
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modifications that tend to mitigate the roundoff errors that plague 
GQR analysis of nonsymmetric distributions. 

II. THEORY AND PROPERTIES OF GQR 
2. 1 Classical use of GQR 

The gqr has traditionally been used as a numerical integration 
procedure and is particularly efficient for computing integrals of the 
form 



r 



f{x)w(x)dx 



where the integrand has been factored into a non-negative term w{x) 
and a strongly continuous term f{x). 

The first application of gqr in the communications literature 1 was 
motivated by the work of Golub and Welsch 2 and Sack and Donovan 3 , 
who showed that the non-negative factor w(x) need not even be 
completely known to compute the desired integral. Only a finite 
number of the moments of w(x) are required to find the desired 
integration rule. Benedetto et al. 1 noticed that the problem of error 
probability evaluation in the presence of intersymbol interference (isi) 
could be posed in this form. Other applications of the gqr technique 
can be found in Refs. 4 through 9. 

In this paper, we apply the GQR technique to a larger class of 
problems where f(x) need not be continuous. We begin by reviewing 
a fundamental result in the theory of gqr. 

Theorem: Let w(x) be a non-negative weight function defined on (a, 
b). Then iff{x) has continuous derivatives up to order 2N (see Refs. 
10 through 13), 



-r 



f(x)w(x)dx 



where 



N 

= % Aifiti) + RM& a<i<b 

t=i 

a < U < b i = 1, 2 - - • N, (1) 



f {m (0 



f {2m (x) is the 2Nth derivative off(x) and (2N)\ is 2N factorial. The 
nodes {%} are the distinct real roots of the unique Nth degree 
polynomial 
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Pn{x) = Aa/II (x- tt), k N >0. (3) 

The polynomials p n (x) are orthonormal with respect to w{x), i.e., 



w(x)p m (x)p n (x)dx = | Q 



m = n 
m^ n. 



The strictly positive weights (or Christoffel numbers) are in turn 
given by 

Ai = Z T^ TTT^TTT i = l > 2 ' "*■ (4) 

&N OjV+1 (£,)/? N (ti) 

where 

dp N (t) 



pWi) = 



dt 



The 2N- tuple {y4„ £,}£,i is known as the iV-point rule corresponding to 
w(x). 

If f(x) is a polynomial of degree (2N — 1) or less, the remainder 
Rn(£) equals zero and the gqr is exact. This affords the maximum 
degree of precision (i.e., the maximum degree polynomial that is 
integrable with no error for an iV-point rule) possible with a quadrature 
formula of the form of (l). 10 " 12 When the remainder is not zero, it can 
be bounded in magnitude to obtain upper and lower bounds on J. The 
bounds obtained in Ref. 1 for the isi and Gaussian noise problem are 
often loose though, and convergence of the iV-term summation in (2) 
is usually much faster than might be inferred from bounds on Rh{£). 

2.2 Methods for computing GQR 

Several algorithms are known for efficiently computing the rule for 
an arbitrary weight function w(x). Extremely useful procedures have 
been discovered by Golub and Welsh, 2 Sack and Donovan, 3 and 
Gautchi. 14 The outstanding merit of these techniques is that the Ap- 
point rule corresponding to a given w(x) can be computed from the 
moments 






x'w{x)dx. (5) 



Because explicit knowledge of the weight function w(x) is not required, 
the gqr procedure is a powerful tool for the analysis of communications 
systems. 

Details of an algorithm for computing gqr are given in the appendix. 
Our algorithm is a modification of Gautchi's procedure, M which tends 
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to reduce computer roundoff errors. The critical stage in the algorithm 
is the Cholesky decomposition of a positive definite matrix of moments. 
The standard Cholesky decomposition used in Refs. 1 and 2 fails when, 
because of limitations of machine accuracy, the matrix is no longer 
positive definite due to roundoff errors. Improved accuracy is obtained 
by using an alternate method of performing the Cholesky decomposi- 
tion* that avoids taking a square root at each step in the algorithm. 15 
Combining the alternate Cholesky decomposition with the modified 
moment algorithm of Gautchi 14 yields an extremely stable method for 
obtaining gqr. Further discussion of techniques to mitigate computer 
roundoff errors is found in the appendix. 

2.3 Computing the distribution of a random variable via GQR 

In Ref. 1 gqrs are used to obtain the exact probability of error for 
digital transmission in the presence of isi and Gaussian noise. The 
problem was reduced, via the gqr approach, to computing the mo- 
ments of the isi and letting f(x) in (1) be the probability of error 
caused by Gaussian noise conditioned on the isi. The isi moments can 
be computed via Prabhu's method 16 when the data symbols are inde- 
pendent. For a large class of correlated data, the moments can be 
efficiently computed via the modified Cariolaro-Pupolin algorithm. 17 ' 18 
Both of the above procedures are easily implemented and have a 
complexity that grows only linearly with pulse duration. 

While there have been numerous applications of GQR to problems in 
the literature, all those known to us have had the restriction that the 
function f(x) has continuous derivatives up to order 2N. Presumably, 
this is because of the desire for strict bounds on the error term in (2). 
If we are willing to forego the analytical error term and consequently 
accept an empirical convergence of (1), we can apply the gqr technique 
to a larger class of problems with excellent results. 

The following theorem shows that no continuity requirements need 
be imposed on f(x). 

Theorem: (see Ref. 19) If W(x) is a fixed, nondecreasing function with 
infinitely many points of increase and the Riemann-Stieltjes integral 

-b 
f(x)dW(x) 



J a 



exists, then 



f 

J a 



N 

f(x)dW(x)= lim ^Aifiti), (6) 



* Applying the alternate Cholesky decomposition to the gqr problem was suggested 
by L. Kaufman. Subsequently, the same approach was found to have been independently 
proposed in Ref. 4. 
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where {Ai, t,)$Li is the gqr corresponding to the moments 
•t, 

l dW(x) i = 0, 1, 2, . ■ • 2N. 



tn = x\ 

J (I 



The function f{x) is arbitrary as long as the integral in (6) exists. 

Because the cdf of a random variable is a nondecreasing function, 
we write the statistical expectation of the function f(x) as 

E[f(x)]= I f(x)dW(x), (7) 



where W(x) is a probability measure with infinitely many points of 
rise. Choosing /(jc) to be the indicator function 

f(x) = $ a (x) 

1 x < a 



(8) 
x > a, 

we obtain the distribution function of the random variable as 
$ a {x)dW(x) = lim X At 



f 



where 



and 



= UmW N (a), (9) 



W N (a) = £A l (10) 



S%= {i\ti^a} 

is the set of indices for which t, < a. 

Since the rule can be obtained from the {/x,-}f^o, we have a means of 
constructing an approximation to the cdf of a random variable from 
its moments. In the limit as N approaches infinity, eq. (9) is exact at 
each point a. 

This leads us to propose the following estimator 

W(x) = W N (x) = X Ai. (11) 

This estimator gives a staircase approximation to the true cumulative 
distribution that becomes increasingly fine as N increases. Equiva- 
lently, each (At, U) can be considered a point mass of a discrete 
approximation to the true probability density function. 
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While Szego's theorem proves the asymptotic convergence of the 
estimator to the true cdf when W{x) has an infinite number of points 
of rise, a different result holds for a discrete distribution with a finite 
number of points of increase. 

Theorem: If W(x) is a fixed, nondeereasing function with M < o° 
points of increase, then 

$ a (x)dW(x) = lim W w (a). (12) 



r 



Proof: An alternate formulation of Gaussian Quadrature 14 is as the 
purely algebraic solution to 

fij = 2 MttV j = 0, 1, ■ 2M. (13) 

Now we assume the unknown discrete PDF is of the form 

M 

w(x)= %AMx-ti). 

1=1 
The moments of this random variable are given by 

M 

which is identical to (13) for N = M. 

Finally, we consider the behavior of the gqr for N larger than the 
number of points of increase M. The result is that the algorithm breaks 
down entirely. This is because a discrete distribution that takes on 
exactly M values is completely characterized by its first 2M moments 
and the addition of redundant moments to the problem causes the 
procedure to fail when the Hankel matrix of moments [eq. (21)] 
becomes nonsingular. 

2.4 A modified GQR estimator 

The following is a modification of the estimator Wa?(«) that has been 
found to be more accurate in many applications. Instead of assuming 
that the approximation PDF is composed of point masses, we assume 
that each area of mass Ai is more accurately modeled by a narrow, 
even symmetric, distribution centered around the point U. Thus, we 
propose the smoothed estimator W${<x) which, when evaluated at a 
node, equals 

W%(td = W N (ti)-~. (14) 

Between nodes, W%{a) is given by any "smooth" interpolation routine. 
A simple linear interpolation was found to be sufficient in the examples 
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Fig. 1 — Convergence of gqr estimator for Gaussian pdf. 



that follow. This estimator does not have the jump discontinuities of 
the estimator WnM and is intuitively more satisfying because it fits 
a smoother distribution to W(x). Wjv(a) has been found to give more 
accurate results when applied to known continuous distributions and 
to discrete distributions when M is much greater than N. 

III. APPLICATION TO ARBITRARY DISTRIBUTIONS 
3. 1 Known continuous random variable case 

To show the convergence properties of the gqr technique, we 
illustrate the behavior of Wn(oc) and Wt/(a) with some examples. We 
begin with the Gaussian distribution. Assuming a zero mean, unit 
variance random variable X, we compute the gqr estimators for 
various values of N in Fig. 1. Reasonably accurate results were obtained 
at the 10" 6 point for N > 10. This empirical rate of convergence is also 
typical of distributions that have near- Gaussian statistics. The GQR 
algorithm, using the Cholesky decomposition described in the appen- 
dix, returned accurate results for all N < 60, where N = 60 was the 
dimensionality limit in the computer program. 

In general, the gqr algorithm performs well for zero mean, symmet- 
ric distributions. To illustrate the problems that can occur with non- 
symmetric distributions, consider the lognormal distribution related to 
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the Gaussian distribution by Y = e x . Straightforward computation of 
the moments yields 

^ = exp(A 2 /2). (15) 

Using these moments in the gqr algorithm, the algorithm breaks down 
at N = 17 because of roundoff errors in computation. This problem is 
solved by a transformation that symmetrizes the distribution. We then 
compute the gqr corresponding to the symmetrized distribution and 
take the inverse transform to obtain the original distribution. 
For the lognormal distribution, we form a new PDF 

w e (y) = \/2\w{y) + w(-y)l (16) 

which corresponds to the even part of w(y). The moments of w e ( y) 
are obtained by setting the odd moments of w{y) equal to zero. The 
symmetric moments are then used in the gqr algorithm to obtain 
W e (y), which is easily transformed back to the desired cdf via the 
relation 

Experience with these procedures suggests that it is well worth the 
effort to transform distributions that are not symmetric and even (see 
Fig. 2). The modified moments produce the same robust accuracy seen 
with the Gaussian distribution above. 

As another example, consider a uniform distribution defined on the 
interval (-1, 1). The convergence of the gqr estimator Wn(<x) is shown 
in Fig. 3. Since the distribution has only finite support, by eq. (1) we 
know that all the nodes will lie in the interval (-1, 1). In the limit as 
N _» oo, the nodes will become more densely packed in this interval 
and 

lim{max|fc|} = 1. (18) 

Thus, the gqr algorithm can be used to find the maximum value that 
a random variable attains, i.e., the largest node t max . This can be used, 
for example, to find the maximum eye degradation in a digital regen- 
erator caused by correlated intersymbol interference. 

All the examples so far have been trivial applications since we knew 
the real distributions a priori. A more interesting application is deter- 
mining the distribution of the sum of K lognormal random variables. 
This problem has a long history and no closed form solution is known. 
This PDF is related to the distribution of crosstalk power in paired 
cable transmission systems and also results from transmission over 
certain types of fading channels. Utilizing the gqr technique, we can 
find the desired distribution if we can compute the necessary moments. 
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Fig. 2— Convergence of gqr estimator for lognormal PDF. 
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Fig. 3 — Convergence of gqr estimator for uniform PDF. 
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Fig. 4— Sum of K lognormal (K = 1, 2, 4, 8). 



Assuming that the lognormal random variables are independent, and 
following Prabhu, 16 we find that the moments of 

V K = J Yi 



are given by the recurrence relation 

E\_(V k Y] = i Q E[(V k -d l ]fH-i, 



(19) 



where {ju.,} ?*i are the moments of the independent, identically distrib- 
uted lognormal random variables. Figure 4 shows the resulting distri- 
butions for K = 2, 4, 8, and 16. As we mentioned above, the distribution 
was symmetrized and inverse transformed to reduce the effects of 
roundoff errors. This technique can be applied to any number of 
arbitrary distributions for which the required moments can be com- 
puted. 

3.2 Known discrete random variable case 

In this section, we apply the gqr estimator W%(a) to discrete 
distributions. First we consider the case of a mixed distribution com- 
posed of a Gaussian distribution plus discrete components. The weights 
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Fig. 5 — Mixed Gaussian distribution, 

and nodes of the estimator are shown in Fig. 5, where the second 
moment of the discrete part equals ten times that of the Gaussian 
component. As we can readily discern, the gqr procedure is useful in 
identifying the discrete components of a PDF. 

As the final example of a known distribution, we consider the sum 
of nine equally spaced delta functions 

w(x) = 1/9 £ S(x - Xi) xt = -5 + i i - 1, 2, ... 9. (20) 

The convergence of the gqr estimator W%(x) is shown in Fig. 6, where 
the N = 9 estimator is exact since the distribution is uniquely defined 
by the first 2N = 18 moments. For N > 9, the algorithm breaks down. 

IV. SUMMARY 

An estimator based on gqr has been proposed, which converges 
rapidly to the cdf of a random variable and requires only knowledge 
of the moments of the random variable in question. The technique is 
generally applicable to a large class of communications problems and 
provides a practical solution to many analytically intractable problems. 
The technique works equally well for discrete and continuous distri- 
butions and assumes no a priori knowledge of the distribution. 
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APPENDIX 

Details of the GQR Algorithm 

In this appendix, we outline the algorithm used to compute Gaussian 
Quadrature Rules. The procedure combines Gautchi's modified mo- 
ment technique 14 with the Cholesky decomposition suggested by Mar- 
tin et al. 15 The resulting algorithm has been implemented using double 
precision arithmetic and has proven stable and robust. 

To compute the 2N unknowns {Ai}fL) and {£,}<^i> we first form the 
matrix of modified moments 

"mi.i mi,2- ■ ■ nn,N+i 
m-2,] 



M = 

TTl2,N'+\ • ■ • mN+l,N+\ 

where my is given by the inner product 

Wly = (Tj-l, Tj-l) 



(21) 



•ft 

i Ti-i(x)Tj-i(x)dW(x) ij = 1, 2 ■ - ■ N + 1 (22) 
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r 



and {Ti) £o are the first N + 1 members of an arbitrary set of 
polynomials satisfying the recurrence relation 

xTj(x) = ajT J+1 (x) + bjTM) + CjTj-^x) j « 0, 1, 2, • - • N 
T_,(x) = 0, aj¥=0. (23) 

The orthogonal Tchebycheff polynomials determined by (23) consti- 
tute a convenient choice, with 

ao = 1 

aj = Cj=Vi j = 1, 2 ■ - ■ 

bj=0 >=1,2-... (24) 

The modified moments m,j in (23) are simply linear combinations of 

the moments 

rb 



f 

•J a 



and can be simplified for the case of the Tchebycheff polynomials by 
using the relation 

Ti(x)Tj(x) = l HT i+J (x) + T H {x)) i >>. 

Thus, if we define 

T k (x)dW(x), 



(25) 



n 



-\ 



then 



rtiij = l h{vi+j-2 + v,-j} i^j. 



(26) 



It is not necessary in theory for the T k {x) to be orthogonal. The 
formulation by Golub and Welsch 2 used the unmodified moments 
corresponding to T*(x) = x h and, hence, % = 1, bj = 0, and c, = for all 
/. As Gautchi shows, 14 the use of modified moments results in less 
sensitivity to computer roundoff errors. 

We next form the tridiagonal matrix 



J = 



'ai ft 

jSl «2 /?2 
ft - 












(27) 



where 



aj = bj + 



r iJ+i 






#, a, 



7-1,2 •■■ N- 1. 
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The r y are found from the relation 

M = R T R. (28) 

The matrix R is an upper triangular matrix and theoretically is 
positive definite if AT is positive definite. In practice, however, M can 
be ill-conditioned and finite precision arithmetic will cause the matrix 
to appear singular. 

The elements of R are related to the moment matrix M by the 
relations 

\l/2 



r « = I nm - £ r\i 1 

r v = 1 rriij - £ r ki r kj J / r u i < j 



i,j= 1,2 ... JV. (29) 

In practice, the computation of R from (29) will fail at relatively small 
values of N when the square root of a negative number is attempted. 
A refined Cholesky decomposition 15 overcomes this problem by only 
requiring square roots to be computed at the end of the decomposition 
and not at each step as in (29). If we define R* by the relation 

R = R*diag(r„), 
then R* will be a unit upper triangular matrix and 

M=R T R 

= R* T diag(r%)R* 

= R* T DR*, (30) 

where D is a positive diagonal matrix. Then, defining the auxiliary 
quantities 

ml = rf/dj, (31) 

the following solution is obtained 

y-i 
mfj = m y - £ m%r% > = 1, 2 • • ■ i - 1 
jt-i 

di — m„ - X m&rft. (32) 

k-i 

The advantage of the alternate decomposition is that square roots 
are not required until the final step, when the positive diagonal matrix 
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Table I — Comparison of three implementations 







Alter- 


Alternate 




Stand- 


nate 


Cholesky 




ard Cho- 


Cho- 


Modified 




lesky 


lesky 


Moments 


Gaussian Random Variable 


17 


60 


60 


Uniform Random Variable 


13 


40 


38 


Lognormal Random Variable 


— 


17 


14 


Symmetrized, Lognormal Random Variable 


— 


60 


60 



D in (30) is factored. Along with the modified moment procedure, the 
alternate Cholesky decomposition yields accurate results even for large 
values of N. 

Several implementations of the gqr algorithm have been examined 
to elucidate the features that contribute to the reduction of computer 
roundoff errors. These include: 

(i) Standard Cholesky 
(ii) Alternate Cholesky 

(Hi) Alternate Cholesky with modified moments 
(iv) All of the above using symmetrized moments. 
Each approach was evaluated in double precision arithmetic. 

The value of iV at which the Cholesky decomposition fails was 
chosen as the measure of robustness for a variety of input probability 
density functions. Some of these results are tabulated in Table I. The 
standard Cholesky consistently had the poorest performance for all of 
the distributions considered. For symmetric distributions, the alternate 
Cholesky scheme provided a significant reduction of computer error. 
For the Gaussian distribution, the procedure was accurate for all 
JV < 60, where 60 was the dimensionality limit imposed on the com- 
puter routine by storage requirements. The addition of the modified 
moment approach resulted in virtually no improvement relative to the 
alternate Cholesky implementation alone. None of the first three 
approaches proved satisfactory for nonsymmetric distributions (e.g., 
lognormal). The solution to this obstacle for one-sided distributions is 
to symmetrize the distribution according to (16), find the gqr estimate 
for the symmetrized distribution, and then obtain the desired distri- 
bution using (17). As we see in Table I, this renders the lognormal 
estimate as robust as the symmetric Gaussian distribution. 

The final step in obtaining the nodes and weights involves finding 
the eigenvalues and eigenvectors of the matrix J in (27). The eigen- 
vector q } corresponding to the eigenvalue tj is found from the equation 



jQj = tjQj j = 1, 2, 



N. 



(33) 



The eigenvalues {£/}w are the nodes of the gqr and the positive 
weights are given by 
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where 
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FORM 
MATRIX J 
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COMPUTE 

(tiMa,) 



UNDO 
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I 



COMPUTE 
DISTRIBUTION 



Fig. 7 — Flowchart of GQR algorithm. 
Aj = q%po, 

qj = (qy, q 2j • ■ • QnJ ■ 



(34) 



A flowchart of the steps used to compute gqr is shown in Fig. 7. 
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